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Reduced density matrix hybrid approach: Application to electronic energy 
transfer 

Timothy C. Berkelbach/'B Thomas E. Markland,^ and David R. Reichman-^ 

Department of Chemistry, Columbia University, 3000 Broadway, New York, New York 10027, 
USA 

Department of Chemistry, Stanford University, 333 Campus Drive, Stanford, California 94305, 
USA 

Electronic energy transfer in the condensed phase, such as that occurring in photosynthetic complexes, fre- 
quently occurs in regimes where the energy scales of the system and environment are similar. This situation 
provides a challenge to theoretical investigation since most approaches are accurate only when a certain en- 
ergetic parameter is small compared to others in the problem. Here we show that in these difficult regimes, 
the Ehrenfest approach provides a good starting point for a dynamical description of the energy transfer 
process due to its ability to accurately treat coupling to slow environmental modes. To further improve on 
the accuracy of the Ehrenfest approach, we use our reduced density matrix hybrid framework to treat the 
faster environmental modes quantum mechanically, at the level of a perturbative master equation. This com- 
bined approach is shown to provide an efficient and quantitative description of electronic energy transfer in 
a model dimer and the Fenna-Matthews- Olson complex and is used to investigate the effect of environmental 
preparation on the resulting dynamics. 



I. INTRODUCTION 

Recent experimental observations of long-lived elec- 
tronic coherence in photosynthetic complexes and solu- 
tions of conjugated polymers^^^ have challenged the con- 
ventional view that environmental effects rapidly quench 
quantum coherence at room temperature. These experi- 
ments have spurred the theoretical investigation of elec- 
tronic energy transfer (EET), including the development 
of accurate numerical methodologies. Such work has 
sought to explain the origin of the observed coherence 
lifetime, to predict the effects of system parameters, and 
to ultimately understand the role of quantum coherence 
in promoting or inhibiting efficient energy transfer. 

EET systems typically consist of a set of molecular 
chromophores which are electronically coupled to one an- 
other as well as to local environmental phonons. The 
former couplings facilitate exciton delocalization while 
the latter couplings tend to destroy this so-called quan- 
tum coherence. Hence, accurately modeling the inter- 
play between these effects is crucial. A difficulty which 
arises in the modeling of EET is the similarity of energy 
scales describing these competing phenomena. For ex- 
ample, in typical multi-chromophoric systems, the elec- 
tronic couplings, reorganization energies, and character- 
istic environmental frequencies are all on the order of 
10 - 100 cm-^ 

As examples of perturbative methods which may fail in 
this regime, the popular Redfield^^^^^ and Forster— ^ the- 
ories of energy transfer implicitly require nearly Marko- 
vian, non-adiabatic dynamics characterized by environ- 
mental frequencies much larger than either the reorga- 
nization energy or electronic couplings, respectively. It 
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is one of the goals of this work to show that the viola- 
tion of this expectation in the form of atypically small 
environmental frequencies actually suggests a promising 
route towards accurate theories of EET in this intermedi- 
ate coupHng regime. Specifically, quantum-classical the- 
ories are well suited to the problem of EET, where en- 
vironmental fluctuations take place over timescales on 
the order of the electronic motion. In such approaches, 
the environmental degrees of freedom are treated classi- 
cally since quantum effects are expected to be insignif- 
icant for such low-frequency vibrational motion. Most 
relevant to EET, these methods can accurately describe 
non-Markovian effects such as transport mediated by 
non-equilibrium phonon states. Investigations along this 
line have included the application of LSC-IVR,^^ the 
Poisson bracket mapping formalism, iterative linearized 
propagation schemes, and a modified variant of mean- 
field Ehrenfest dynamics. In a similar vein, the re- 
cently developed reduced hierarchy equations,— although 
quantum-mechanically exact when fully converged, are 
numerically simplest for systems at high-temperature or 
with a slow, adiabatic bath. 

In accord with the above discussion, in the present 
work we show that a simple quantum-classical Ehren- 
fest treatment of EET systems yields results in quali- 
tative and even sometimes semi-quantitative agreement 
with existing exact results. In particular, quantitative 
agreement is found for short- time dynamics, as well as 
coherence frequencies and lifetimes. This latter finding 
provides insight into the mystery of long-lived quantum- 
coherence: the intramolecular motions are simply too 
slow to induce effective dephasing. We also demonstrate 
that including the quantum subsystem's back-reaction on 
the classical environment improves the long-time popula- 
tions in contrast compared to a previous Ehrenfest study 
where this effect was not included. However, although 
this improves the accuracy of long-time populations, the 
Ehrenfest method generally still yields incorrect long- 
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time values due to a well-known intrinsic violation of de- 
tailed balance. These long-time properties are crucial for 
the accurate description of EET, where the population 
of target sites provides a simple metric for the overall 
transport efficiency. 

To address this issue, we invoke our recently intro- 
duced reduced density matrix hybrid (RDM-Hybrid) 
methodology^^ which enables one to treat the system and 
high frequency environmental modes quantum mechani- 
cally, while the slow modes are handled by the Ehren- 
fest approach. This methodology allows for efficient 
treatment of quantum and classical environmental modes 
and is able to quantitatively correct the discrepancies of 
Ehrenfest dynamics, yielding excellent agreement with 
the exact results obtained by Ishizaki and Flemin g^^i^^ 
using the reduced hierarchy equations (RHE). As an ex- 
ample of the types of problems one may further probe 
using our physically transparent methodology, we inves- 
tigate the effects of the initial bath preparation on the 
subsequent observation of coherent quantum dynamics. 

The outline of the paper is as follows. In Sec. [Ill we 
present the Frenkel exciton Hamiltonian used for EET 
modeling. We then review the Ehrenfest and RDM- 
Hybrid methodologies in Sec. IHII In Sec. |lVl we de- 
rive a perturbative quantum master equation to use in 
our RDM-Hybrid algorithm for the system and quantum 
environment modes. We present Ehrenfest and RDM- 
Hybrid results in Sec. |Vl including population dynamics 
and rate constants for a simple dimer, as well as popula- 
tion dynamics for the seven-site Fenna-Matthews-Olson 
complex with two types of environmental preparation. 
We conclude in Sec. IVTl 



ground state. The electronic coupling of this excitation 
between sites m and n is given by Jmn and is assumed 
to be static. 

In contrast to the spin-boson Hamiltonian investigated 
in our previous work, which described a two-level system 
coupled to one shared bath,^^ the present Hamiltonian is 
more commonly adopted for molecular energy transfer, 
where molecular vibrations and surrounding environmen- 
tal effects are local and approximately uncorrelated. 

The phonon bath located at each site is completely 
characterized by its spectral density, which we take to be 
of the Debye form. 



In the above, ujc is the bath cutoff frequency and the reor- 
ganization energy, A = 7r~^ duo J {uj) / uo ^ is the energy 
dissipated by the environment after a Franck- Condon 
transition of the excitation from one site to another. In 
the second line and henceforth, we assume that all sites 
have the same spectral density. However we note that the 
methods employed here are in no way limited (or catered) 
to the treatment of independent baths, the functional 
form of the spectral density, or the assumption of iden- 
tical baths. Future work will include the investigation of 
each of these effects on energy transfer dynamics. 



II. MODEL HAMILTONIAN 

As in other theoretical work on EET, we adopt the 
Frenkel exciton Hamiltonian for N chromophores, given 

by 

N N 

H=Y^ \n)En{n\ + ^ \m)J^n{n\ + Hb + H^b (1) 

n=l m^n 

where one assumes independent baths for each site, 

N 

Hb = J2T.[Kn+<nQln]/2, (2) 

n—1 k 

and a bi-linear system-bath coupling, 

N 

Hsb = ^ \n){n\^Ck,nQk,n' (3) 

n=l k 

Such a Hamiltonian physically describes the single- 
excitation subspace of the complex's total Hilbert space, 
where En denotes the energy of the total system when 
the nth chromophore is excited and all others are in their 



III. REVIEW OF EHRENFEST AND RDM-HYBRID 
METHODS 

In this section, we briefly review the pertinent de- 
tails of the Ehrenfest and reduced density matrix hy- 
brid (RDM-Hybrid) quantum dynamics methodologies 
needed to treat to the Hamiltonian in Eqs. ([I])- ([3]). For 
full details of the method, we refer the reader to our pre- 
vious paper 



A. Ehrenfest method 

In the Ehrenfest method, one assumes separability of 
system and bath variables, yielding the product density 
matrix p{t) ^ ps{t)pb{t). Inserting this ansatz into the 
Liouville equation and tracing out the system or bath 
variables yields the time-dependent self-consistent field 
coupled equations of motion. A classical treatment of the 
bath density operator then yields the Ehrenfest equations 
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of motion, 

dp.s{t) 
dt 



-i[H,,p,{t)] 

N 



\n){n\^Ck,nQk,n{t),Ps{t) 



.n=l 



-w| „Qfc,„ - Ck,nT£s{\n){n\ps{t)} 



(5) 

(6) 
(7) 



dQk,n 

dt 

dPk,n 

dt 

where square brackets denote the commutator. The fluc- 
tuating bath coordinates yield a time-dependent bias for 
the system Hamiltonian, and in-turn the system pop- 
ulations TTs{\n){n\ps{t)} = Pn{t) yield a driving force 
for the bath coordinates. The system populations Pn{t) 
are averaged over an ensemble of such coupled trajecto- 
ries with the harmonic bath initial conditions sampled 
from the classical Boltzmann or the quantum Wigner 
distribution.^^ 



B. RDM-Hybrid method 

The RDM-Hybrid method differs from the Ehrenfest 
approach by partitioning the bath degrees of freedom into 
'core' and 'reservoir' modes. In terms of the spectral 
density, we have 

JcoreM= (8) 
Jres(^) = J{UJ)S{U,UJ*), (9) 

where S{uj^uj^) is a switching function, taken here to be 



(10) 



which switches smoothly from 1 to as goes from to 
cj*. The switching frequency, cj* is taken to be a charac- 
teristic timescale of the electronic system^ 

Positing the separation of 'system-core' and 'reservoir' 
density operators, p{t) ~ psc{t)Pr{t), one finds a Liouville 
equation for the density matrix of the system and core 
modes, Psc{t), given by 



dpscjt) 
dt 



-i [Hsc{t),psc(t)] , 



(11) 



where the time-dependent Hamiltonian, Hsc(t)^ is a mod- 
ified system-core Hamiltonian, 

N N 

Hsc{t) = \n)En{t){n\ + ^ \m)Jmn{n\ 

n=l m^n 



N 



n—1 /cGcore 



\n){n\ck^nQk,', 



(12) 



and the time-dependent bias arises from the coupling to 
the classical reservoir modes, 



En{t) = En+ J2 (^kQkit). 



(13) 



feGres 



A solution of the above Liouville equation for the total 
system-core density matrix is of course intractable, but 
its reduced density matrix averaged over the core degrees 
of freedom can be calculated by a variety of existing ap- 
proximate and exact methods. The diagonal elements 
(populations) of this reduced density matrix, 

TisTTc{\n){n\psc{t)] = TTs{\n){n\ps{t)} = Pn{t), (14) 

in turn yield a driving force in the classical reservoir equa- 
tions of motion. 



dQk.n 

dt 
dPk.n 
dt 



•^k,n 



Qk,n (^k,nPn{t) ' 



(15) 
(16) 



As in the Ehrenfest method, the final system popula- 
tions are calculated as an average over trajectories. In 
this work, the initial conditions of the reservoir modes 
are sampled from the Wigner distribution. To evolve 
the reduced system-core density matrix we perform an 
approximate calculation using a perturbative quantum 
master equation. This execution of the RDM-Hybrid ap- 
proach yields an efficient methodology not much more 
expensive than a typical master equation calculation but 
with far superior accuracy 

As shown in our previous paper, the RDM-Hybrid 
methodology naturally interpolates between the regimes 
of validity of its composite methods and furthermore 
works well even in regimes where neither method alone is 
accurate. By employing the Ehrenfest method, accurate 
for nearly adiabatic dynamics, and a non-adiabatic quan- 
tum master equation, the RDM-Hybrid approach can ac- 
curately treat the entirety of parameter space using a 
single dynamical scheme. The perturbative master equa- 
tion used in this work to treat the driven system-core 
dynamics is derived in the following section. 



IV. MASTER EQUATION FOR SYSTEM AND CORE 

Following the success of our previous work,^- 
which employed the noninteracting blip approximation 
(NIBA)^i2i for the system -core reduced dynamics, we 
here derive a multi-site generalization, perturbative to 
second order in the electronic couplings, Jmn, which 
we shall continue to refer to as 'NIB A,' for simplicity. 
A similar master equation, the noninteracting cluster 
approximation,^^ has been derived by different means. 

Following Hu and Mukamel^^ (see also Golosov and 
Reichman^^), we define a Liouville-space projection op- 



4 



erator, 



N 



J2\np,)){{n\. 



(17) 



Here |n)) = |n)l5(n|, I5 is the identity operator in the 
bath degrees of freedom, \npi))) = |n)p5(n|, and the in- 
ner product is given by {{A\B)) = TTsTn{A^ B). In this 
notation, the observables of interest, i.e. the site popula- 
tions, are written as Pn{t) = ((n|p(t))), where p{t) is the 
total system-bath density matrix. 

Employing the usual projection operator 
formalism^— one may derive the exact set of 
equations 



Pn{t) 



Si 



dTKnm(t,T)Pm{T) 



where the kernels are given by 

Knmit,T) = {{n\Cv{t)U{t,T)Cv{T)\mpb)). 



(18) 



(19) 



In the above expression, Cv{t) • • • = [V{t)^ . . .] is the 
Liouvillian in the interaction picture, i.e. V{t) = 
ex.p{iHot)V ex.-p{—iHot). The propagator is given by 



U (t, r) = exp^ 



dt'QCvit') 



(20) 



with exp^[. . . ] denoting the usual time-ordered exponen- 
tial and Q = 1 — P is the complementary projection op- 
erator. The perturbation here is ^ = J^^-^^ |^)^mn(^| 
with the unperturbed Hamiltonian given simply by Hq = 
H-V. 

While the above formalism is exact, the propagator 
U{t,r) in the complementary subspace is intractable. 
However the propagator may be expanded perturba- 
tively, and to lowest non-trivial order in the electronic 
couplings (obtained by setting U{t,r) = 1), one obtains 
the second-order kernels, 

KnM^^ ^) = 2^nm ^xp [-Q2{t - t)] 

X cos|cn(^,r) -Cm(^,r) -^Qi{t-r) 

^^^^^\Qi(t)-Qi(r)]l (21) 



with Knn{t,r) = -Y.m^n^rnn{t,r). In the above, we 
have introduced the site-dependent accumulated phase. 



Cn(t,r) 



dt'Enit') 



= / dt' 



^ ^ Ck,nQ k,n(t ) 
/cGres 



(22) 



and the bath correlation function, Q{t) = Q2{t) + ^Qi(t), 
is given by 



Qit) 



2 r 

TT Jo 



^^ Jcore(^) {(,oth(/3w/2) [1 - COs(wt)] 



+ i sin(wt)} . 



(23) 



The factor of two appearing here differs from the usual 
NIBA expressions due to the assumption of independent 
baths. Lastly, the shift parameters Sn arise from the 
initial bath density matrix. 



Pb 



n ^ exp ( -P hk,n 

n \ k J 



where 



hk,n = - Pk,n + ^fc,n {Qk,n - ^nCk.nf' 



(24) 



(25) 



We now demonstrate that our generalized master equa- 
tion, Eq. (dH), naturally reduces to the hopping rate 
equation predicted by Forster theory in its regime of 
validity. Specifically, we consider the strongly non- 
adiabatic regime, such that none of the bath modes are 
treated classically, i.e. (n{t^r) = En x {t — r). For sim- 
plicity we consider the case Sn = = —I such that the 
final term in Eq. (|2T]) vanishes. 

When the bath dynamics are much faster than those of 
the electronic subsystem, the memory kernel in Eq. (p!8|) 
decays very quickly, such that the Markovian approxima- 
tion can be made 



N 



(26) 



m—l 



where the rate constant knm is given by knm = 
dtKnm{t)' Introducing the well-known line-shape 
function g{t) = Q{t)/2 — iXt, with A the reorganization 
energy, we can write 



f 

Jo 



dt 



dt'5{t - 1') 



X e 



-i{Er,+\)t-g{t) i{Em-\)t' -g{t') 



(27) 



Using the Fourier transform of the Dirac delta function, 
5(t — t') = (27r)~-^ do; exp[ia;(t — t')]^ we have 



t2 ro 

27r J_, 



dujAn{uj)Fraiuj) 



(28) 



with the absorption and fluorescence spectra given by 

/oo 
^^^iu.t^-i(E^ + \)t-g(t)^ (29) 
-OO 

/oo 
^^gi^*g-i(B„-A)t-g-(0, (30) 
-OO 
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Note we have made use of the fact that A{uj) and F{uj) 
are real, as can be exphcitly checked by symmetry. Equa- 
tion ([28|) is readily recognized as the celebrated Forster 
rate. 



V. RESULTS 

A. Model EET dimer 

We begin our application of the Ehrenfest and RDM- 
Hybrid methodologies with the simplest model of exci- 
tation energy transfer, a dimer. For the sake of com- 
parison, we employ the model Hamiltonian, Eqs. ([I])-(|3j), 
with a Debye spectral density, Eq. (H]), and the parame- 
ters of Ishizaki and Fleming,— for which Redfield theory 
has already been shown to fail badly. In the results to 
follow, we compare to existing results obtained via the 
numerically exact reduced hierarchy equations (RHE)^ 
based on work originally done by Kubo et a L^^i^^ For 
the physically motivated reasons explained in our previ- 
ous paper, we take the the splitting frequency to be 
equal to the Rabi frequency of the electronic subsystem. 



(a) ?i = J / 50 



{b)X = J / b 



uj"^ = ujR = y {El - E2) 



In all of our results we discretized the bath into / ■ 
modes with frequencies and couplings given by 



ujk = (jOc tan 



L2^(^-l/2) 



2 2 J(wfe) 

Cfc = -Wfc— r 

TT p{LOk) 



(31) 
-- 300 

(32) 
(33) 



for k = 1,2,...,/, which can be shown to reproduce 
the reorganization energy exactly. The coupled system- 
bath equations of motion were solved with a second-order 
Runge-Kutta scheme using a timestep of 0.5 fs. For con- 
sistency, the classical degrees of freedom were sampled 
from the Wigner distribution, though for the relatively 
high temperatures considered here, the results are largely 
insensitive to this choice when compared to those ob- 
tained by purely classical Boltzmann sampling. 



1. Strong electronic coupling 

We begin by considering an EET dimer with relatively 
strong electronic coupling J12 = J21 = J = 100 cm~^ 
and energetic bias Ei — E2 = 100 cm~^, such that the 
splitting frequency, Eq. (|3T]) . is cj* ~ 220 cm~^. The 
simulations are performed at room temperature, T = 
300 K, with a bath cutoff frequency of Uc = 53 cm~^. 
The bath correlation time corresponding to this cutoff 
frequency is Tc = oj~^ = 100 fs. Since uJc/ J ^ 0-5, this 
regime can be characterized as being weakly adiabatic. 



Exact (RHE) 
NIBA 
Ehrenfest 
RDM-Hybrid 
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{c)X = J (d) X = 5J 




200 400 600 800 1000 200 400 600 800 1000 
Time t [fs] Time t [fs] 

FIG. 1. Population of site 1 in an EET dimer system con- 
sidered by Ishizaki and Fleming with Ei — E2 — 100 cm~^, 



J 100 cm"\ 



53 cm ^ 



100 fs), andT = 300 K. 



Each site is coupled to its own bath with a Debye spectral 
density. 



The population of site 1 for these parameters, with the 
initial condition p(0) = |1)(1| exp(— is shown in 
Fig. [1] for a reorganization energy varying over more than 
two order of magnitude. As alluded to in the introduc- 
tion, the Ehrenfest dynamics are seen to be qualitatively 
very good, though for intermediate reorganization ener- 
gies, they show a significant deviation in the long-time 
population. This behavior is to be contrasted with that 
observed using the LSC-IVR method, where the agree- 
ment regularly worsened with increasing reorganization 
energy (see Figs. 1 and 2 of Ref. 15). The Ehrenfest and 
LSC-IVR approaches are almost identical for the simple 
Frenkel Hamiltonian considered here, the only difference 
being that the Ehrenfest approach treats the electronic 
degrees of freedom quantum mechanically whereas LSC- 
IVR treats them by a Meyer- Miller mapping to classi- 
cal oscillator s,^^^^ such that all dynamical degrees of 
freedom are treated on equal footing. The results ob- 
served here indicate that such a consistent treatment does 
not necessarily yield more accurate results, especially in 
regimes of strong reorganization energy. In light of these 
failures of LSC-IVR and Ehrenfest approaches, it is cru- 
cial to notice that our RDM-Hybrid methodology yields 
nearly exact population dynamics for all values of the 
reorganization energy, including an accurate treatment 
of long-time dynamics due to the quantum mechanical 
treatment of high frequency environmental modes. 

Moving deeper into the adiabatic regime, we next con- 
sider the same EET dimer but with a smaller cutoff fre- 
quency, UJc = II cm~^ (longer bath correlation time. 




200 400 600 800 1000 200 400 600 800 1000 
Time t [fs] Time t [fs] 

FIG. 2. The same as in Figure [1] but with ujc = 11 cm~^ 
(cj-i = 500 fs). 



1 10 100 1000 

Reorganization energy X [cm ^] 



FIG. 3. Downhill energy transfer rates for an EET dimer in 
the weakly non-adiabatic regime (cJc/J > 1), with Ei — E2 — 
100 cm"\ J = 20 cm"\ ujc 53 cm"^ (cj^^ 100 fs), and 
T = 300 K. 



be seen in our Figs. [T] and [21 the full Ehrenfest treatment, 
including the back-reaction, always yields more accurate 
long-time populations (less than 1/2). 



Tc = 500 fs), such that uJc/ J ^ O-l- Population dynamics 
are shown in Fig. [2] for the same range of reorganiza- 
tion energies as above. For this smaller cutoff frequency, 
Redfield dynamics have been shown to be inaccurate even 
for the smallest reorganization energy considered. How- 
ever, as discussed in Sec. IH quantum-classical methods 
are highly suitable in this strongly adiabatic regime. In- 
deed, the Ehrenfest results presented here and the LSC- 
IVR results of Ref.dH are nearly exact. Again, the RDM- 
Hybrid results are excellent, correcting the minor discrep- 
ancies seen in the long-time populations of the Ehren- 
fest dynamics. The RDM-Hybrid methodology naturally 
'tunes' itself to the more accurate of its two composite 
methods. For example, in going from Fig.[T]to Fig.[2l low- 
ering the bath cutoff frequency further below the splitting 
frequency results in treating a higher percentage of bath 
modes with Ehrenfest dynamics, the more accurate of 
the two methods in this parameter regime. However, as 
we demonstrated in our earlier work,^^ the RDM- Hybrid 
approach performs better than the sum of its parts and 
can also treat regimes where neither NIBA nor Ehrenfest 
dynamics alone would be suitable. 

A recent work^^ also presented Ehrenfest results simi- 
lar to those shown in Figs. [T] and [2] but the back-reaction 
of the quantum system on the classical one, given at 
the end of Eq. ([7]), was neglected. Thus, the classi- 
cal bath harmonic oscillators were isolated and provided 
only a fluctuating bias in the system equations of motion. 
This approximation is akin to well-known Haken-Strobl- 
Reinecker type of master equation which yields equilib- 
rium populations consistent with infinite temperature, 
i.e. both populations go to 1/2 in a dimer regardless of 
the bias, a deficient behavior observed in Ref. As can 



2. Weak electronic coupling 

Although we argue that many of the interesting dy- 
namical results in the recent EET literature can be as- 
cribed to the system's adiabaticity, we close our study 
of Ishizaki and Fleming's EET dimer by investigating a 
mildly non-adiabatic set of parameters, i.e. a dimer with 
weak electronic coupling, to demonstrate the flexibility 
of our RDM-Hybrid approach and continued success of 
Ehrenfest dynamics. The parameters are the same as 
in the previous section, taking uJc = 53 cm~^, but now 
with J = 20 cm~^, such that that adiabaticity ratio is 
(jOc/ J ~ 3. 

Rather than calculating population dynamics as above, 
we determine the downhill energy transfer rate as a func- 
tion of the reorganization energy, A. Due to the small 
electronic coupling, the population dynamics are gener- 
ally well described by an exponential decay and a simple 
fitting procedure yields the uphill and downhill rate con- 
stants. The results of the Ehrenfest and RDM-Hybrid 
approaches are shown in Fig. [3] and can be seen to be in 
almost perfect agreement with the exact RHE results. 
For completeness, we also present results obtained by 
conventional Redfield theory and the pure NIBA-like 
equations derived in Sec. lIVi Redfield theory is perturba- 
tive in XkT /uol and hence it would be expected to break 
down when the reorganization energy A is on the order 
of 0^1/ kT ^ 10 cm~^, an order-of-magnitude prediction 
which is seen to hold unreasonably well in Fig. [3l Fur- 
thermore, we point out that our NIBA-like dynamical 
theory almost exactly reproduces the Forster result pre- 
sented in Ref. [n|, as one would expect due to their close 
formal relation, discussed at the end of Sec. IIVI 
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FIG. 4. Population dynamics of the FMO complex at T = 
300 K, with Tc = cj^^ = 166 fs and bath initial conditions 
sampled from the Wigner distribution (for both Ehrenfest and 
the present hybrid method). The excitation is initially local- 
ized to site 1 in panels (a) and (b) and to site 6 in panels (c) 
and (d). Ehrenfest [(a),(c)] and the RDM-Hybrid approach 
[(b), (d)] are compared to exact results obtained with the re- 
duced hierarchy equations (filled circles). 
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FIG. 5. Population dynamics of the remaining four BChl sites 
not depicted in Figs. [HJc) and (d), i.e. with site 6 initially 
excited. 



electronically-coupled neighbor. 



1. Comparison with existing results 



B. Fenna-Matthews-Olson Complex 

We now proceed to the investigation of the Fenna- 
Matthews-Olsen (FMO) complex, which has been the 
subject of extensive experimental and theoretical inves- 
tigation, especially pertaining to the origin of long-lived 
quantum coherence. Although an eighth bacteriochloro- 
phyll (BChl) chromophore has recently been identified, 
we consider the FMO model Hamiltonian of only seven 
BChl sites, for which numerically exact results existi^ 
and a variety of other approximate methods have been 
tested. The electronic Hamiltonian is taken from Ref. |38| 
and all sites are assumed to have independent, identical 
baths characterized by a Debye spectral density, Eq. |4l 
with A = 35 cm~^. The bath cutoff frequency and tem- 
perature will be varied throughout our investigation. Our 
initial results in this section will serve to demonstrate 
the accuracy of our RDM-Hybrid approach when treating 
such systems while the second part will use the method 
to study the effect of bath preparation on quantum co- 
herent dynamics. 

When moving from a dimer to a multi-site system, it 
becomes less obvious how to choose a characteristic sys- 
tem frequency for the switching frequency, a;*, required 
by our RDM-Hybrid method. For the FMO system, 
we take the switching frequency equal to the Rabi fre- 
quency of the initially excited site and its most strongly 



We begin by comparing to the numerically exact re- 
sults of Ishizaki and Fleming^^ and using the initial con- 
dition from that work, 

p(0) = |no)(no| eM-PH,)/Z, (34) 

where the initial excitation site is taken to be no = 1 
or 6 due to their proximity to the chlorophyll baseplate. 
This sampling of the bath corresponds to a spectroscopic 
initial condition in which the bath initial conditions are 
not equilibrated to the presence of the system. 

As discussed above, the switching frequency is set 
equal to the Rabi frequency of sites 1 and 2, yielding 
cj* = 210 cm~^, or to the Rabi frequency of sites 5 and 
6, yielding cj* = 220 cm~^, for initial excitations at sites 
1 or 6, respectively. In Fig. SI we consider the param- 
eter set investigated by Ishizaki and Fleming for which 
Ehrenfest dynamics are expected to work best, namely a 
very slow, high temperature bath with Tc = oj~^ = 166 
fs and T = 300 K. Despite the relatively weak reorga- 
nization energy, master equations that are perturbative 
in the system-bath coupling (such as Redfield theory) are 
unable to accurately reproduce these dynamics due to the 
system-bath adiabaticity^^ In contrast, the Ehrenfest ap- 
proach performs very well for short times and accurately 
reproduces the coherence frequency, amplitude, and life- 
time. However, just as in our above study of a two- level 
system, we see that the long-time populations deviate 
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FIG. 6. The same as in Fig. [4] but for the longer bath corre- 



lation time, Tc 



50 fs 



from the exact values, a flaw which is impressively reme- 
died with our RDM-Hybrid approach, yielding excellent 
agreement overall. For example, by comparing Figs. lU^c) 
and (d), we see that the RDM- Hybrid dynamics correctly 
reproduce the population inversion completely missed by 
the Ehrenfest dynamics. 

For clarity, we only show population dynamics for three 
of the seven sites in Fig. IH However, the conclusions 
drawn are entirely unchanged for the four remaining sites 
with smaller populations, as shown in Fig.[5]for the initial 
excitation no = 6. 

By shortening the bath correlation time, we expect 
that the performance of Ehrenfest dynamics should de- 
grade when compared to the exact result, since high fre- 
quency bath modes necessitate a quantum treatment. 
Indeed, this expectation is realized in Fig. [6l for which 
Tc = 0J~^ = 50 fs. Figures [6] (a) and (c) show that Ehren- 
fest dynamics again yield qualitatively accurate coher- 
ence lifetimes but incorrect long-time populations, now 
worsened due to the short bath correlation time. In con- 
trast to Ehrenfest dynamics, the quantum- mechanical 
treatment of high-frequency modes in the RDM-Hybrid 
methodology results in excellent performance, again ex- 
emplified by the population inversion in Fig. [6jd). 

The final set of standard conditions considered here 
again has a short bath correlation time Tc = 50 fs but 
at a reduced temperature, T = 77 K corresponding to 
that at which some of the original experiments observ- 
ing quantum coherence were performed i At this low 
temperature, both approximations used in the present 
RDM-Hybrid implementation (NIBA and Ehrenfest) are 
known to worsen. Nonetheless, the RDM-Hybrid popula- 
tion dynamics shown in Fig. [3 are impressively good, and 
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FIG. 7. Population dynamics of the FMO complex at T = 
77 K, with Tc — uJc^ — 50 fs and bath initial conditions sam- 
pled from the Wigner distribution (for both Ehrenfest and 
the present hybrid method). The excitation is initially local- 
ized to site 1 in panels (a) and (b) and to site 6 in panels (c) 
and (d). Ehrenfest [(a),(c)] and the RDM-Hybrid approach 
[(b), (d)] are compared to exact results obtained with the re- 
duced hierarchy equations (filled circles). 



qualitatively much more accurate than those of Ehren- 
fest dynamics. However, both methodologies again make 
excellent prediction of the coherence frequency and life- 
time, which is the experimentally observed phenomenon 
that has generally garnered the most attention. It is also 
worth comparing again to the LSC-IVR work of Ref. UH, 
which presented population dynamics for the FMO com- 
plex in the present parameter regime (tc = 50 fs, T = 77 
K). The LSC-IVR calculation severely underestimate the 
coherence lifetime, yields some negative populations at 
long times, and requires about 5 x 10^ trajectories to 
achieve convergence - almost two order of magnitude 
greater than required by Ehrenfest of RDM-Hybrid ap- 
proaches. 



2. Effects of bath preparation on coherent transport 

The initial condition considered up to this point, 
Eq. (|34|) . is sometimes referred to as a 'spectroscopic' 
preparation, as it corresponds to the physically correct 
initial condition following a rapid excitation from the 
ground state in accordance with the Franck-Condon prin- 
ciple. Therefore, this initial condition is likely to most 
closely resemble the preparation realized in recent in- 
frared spectroscopy experiments. However in biological 
functioning, it is unlikely that this initial condition is 
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physically correct for the dynamics of the FMO complex, 
for example. Recall that the excitation in the FMO com- 
plex is transferred in from the chlorophyll's baseplate. 
This long-range transfer is likely a non-adiabatic process 
and thus requires a fluctuation in the bath coordinates of 
the FMO acceptor site leading to an excited state geom- 
etry, as in the traditional Marcus picture. In light of this 
discussion it is clearly worth investigating to what extent 
oscillatory population dynamics are modified by a non- 
spectroscopic initial condition and whether the observed 
long-lasting quantum beating - in both experiments and 
simulations - is perhaps a product of unphysical spectro- 
scopic initial conditions. 

Despite the recent interest surrounding quantum co- 
herence in energy transfer molecules and materials, very 
similar investigation of electronic and vibrational coher- 
ence in electron transfer reactions began almost twenty 
years ago (see e.g. Refs. I4QH43I) . Most related to our 
present investigation is the theoretical work of Lucke et 
al.^ who investigated the effects of initial bath prepara- 
tion on the possible observance of electronic coherence. 
The principal conclusion of their work was that absence 
of oscillatory population dynamics does not necessarily 
imply dephasing-induced decoherence and that the way 
in which the bath is prepared can suppress or enhance 
oscillations. 

The RDM-Hybrid approach can very naturally treat 
arbitrary bath initial conditions and here we consider 
two different choices for the initial bath density matrix, 
defined in Eqs. (|M|) - (|25]) . What we will term an 'un- 
shifted' initial condition has 5n = for all sites n (spec- 
troscopic preparation), whereas a 'shifted' initial condi- 
tion has 6no = —I where no is the initially occupied sys- 
tem site and Sn = for all other sites n. Note that 
within the RDM-Hybrid scheme, these shift parameters 
manifest in the sampling of the classical reservoir coor- 
dinates as well as in the master equation memory kernel, 
Eq. 

To best exemplify this phenomenon, we will consider 
a parameter regime for which such sensitivity to initial 
conditions is expected to most strongly manifest, namely 
for an adiabatic, low temperature bath. In particular, 
we consider T = 77 K and Tc = oj~^ = 166 fs with the 
reorganization energy to be varied. In Fig.[8fa), we show 
the effect of shifted initial conditions calculated with the 
RDM-Hybrid approach for the standard reorganization 
energy A = 35 cm~^ and the excitation initially local- 
ized to the first BChl site. Clearly, the shifted initial 
condition yields population dynamics with a reduced os- 
cillation amplitude, by about a factor of two. 

Physically, the above effect occurs because the shifted 
initial condition brings the total system into a near- 
eigenstate, in essence trapping the particle in a quasi- 
stationary state with no observable coherent oscillations. 
On the contrary, with an unshifted initial condition, the 
total system is in a superposition of eigenstates, and co- 
herent dynamics will be observed. 

Further exploring this effect, we increase the reorgani- 




FIG. 8. Effects of initial bath preparation on the FMO pop- 
ulation dynamics at T = 77 K with Tc = cj~^ = 166 fs. 
Solid lines depict the unshifted (spectroscopic) initial condi- 
tion, whereas dashed lines depict the shifted (solvated) initial 
condition, for a range of reorganization energies, A (a)-(c). 
Also shown is a schematic diagram (d) of the lower energy 
adiabatic potential energy surface as a function of the gen- 
eralized reaction coordinates, Qi and Q2, to which sites 1 
and 2 are coupled. Solid line arrows correspond to an un- 
shifted initial condition and subsequent excitation dynamics, 
whereas the dotted line arrow corresponds to an initial condi- 
tion shifted to the minimum of site 1, trapping the excitation. 



zation energy in Figs.[8l^b) and (c) to A = 100 cm~^ and 
350 cm~^, respectively. For A = 100 cm~^, the shifted 
initial condition can be seen to almost completely sup- 
press all population oscillations despite predicting almost 
identical population relaxation. Increasing further to the 
very strong A = 350 cm~^, we show that for the shifted 
initial condition, the bath completely traps the excitation 
at the initial site, whereas for the unshifted initial condi- 
tion, the bath can only trap the particle after short-time 
coherence allows significant population transfer. This sit- 
uation is depicted schematically in Fig. [Sfd), where the 
large, solid arrow marks the unshifted initial condition 
and the smaller solid arrows show the possible relaxation 
into neighboring diabatic wells. Relaxation into the dia- 
batic well of site 2 corresponds to population transfer. On 
the other hand, the large, dashed arrow marks an initial 
condition shifted to the minimum of site 1, effectively 
trapping the excitation and preventing any population 
relaxation. These two contrasting population relaxation 
behaviors can be seen clearly in the dynamics of Fig.[8fc). 
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VI. CONCLUSIONS 

To summarize, we have argued that the low frequency 
environmental motions present in many recent models 
of excitation energy transfer (EET) marks a significant 
deviation from the validity regimes of popular existing 
methods, such as the Redfield and Forster theories. In 
light of this observation, we have shown that quantum- 
classical approaches provide a simple and intuitive route 
towards more accurate modeling of EET in intermedi- 
ate coupling regimes. Although the typical Ehrenfest 
method, which treats the electronic subsystem quantum 
mechanically and the environment classically, yields rea- 
sonable results when compared to exact calculations for 
a wide range of EET Hamiltonians considered in the lit- 
erature, its long-time populations suffer from well-known 
unrestricted energy flow, often yielding equilibrium pop- 
ulations corresponding to an infinite temperature result, 
i.e. equal population of all sites. 

To alleviate this problem, we have employed our re- 
cently developed RDM-Hybrid algorithm, which ex- 
tends the usual Ehrenfest method by including high- 
frequency environmental modes into a quantum "core." 
Importantly, for pure system properties, one need only 
calculate the reduced density matrix averaged over the 
quantum core, which can be done approximately but 
accurately for the high-frequency modes included. In 
turn, the remaining "reservoir" modes are treated clas- 
sically and the usual mean-field coupling exists between 
the system-core and reservoir degrees of freedom. Such 
an approach yields excellent results and can be applied 
without modification to a great variety of system-bath 
Hamiltonians, extending well beyond the domain of EET 
parameters. 

In addition to the very favorable comparison with ex- 
isting results, the RDM-Hybrid method was also em- 
ployed for a novel investigation of initial bath prepara- 
tion and its effects on subsequent population dynamics. 
In particular, we showed that the experimentally relevant 
spectroscopic initial conditions often employed in calcu- 
lations may be partly responsible for the unexpectedly 
long-lived quantum coherence. The degree to which this 
behavior manifests in biological functioning will be sen- 
sitive to the way in which excitations enter FMO and 
related complexes, a topic which has received almost no 
attention in the literature but is surely deserving of fur- 
ther investigation. 

Lastly, we point out that whereas most existing ex- 
act methodologies, including influence functional-based 
path-integral methods and the reduced hierarchy equa- 
tions, are reliant on harmonic bath degrees of freedom, 
the Ehrenfest approximation is applicable for systems 
with generically anharmonic degrees of freedom. Thus, 
the work performed here marks an important step to- 
wards the treatment of EET in realistic molecular sys- 
tems and could easily incorporate low-frequency anhar- 
monic modes. High-frequency anharmonicities could also 
be treated, but would require numerical evaluation of 



the memory kernel using e.g. semi-classical methods4^i^ 
Such work is currently in progress. 
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